library(data.table)
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
data.table 1.15.4 using 28 threads (see ?getDTthreads).  Latest news: r-datatable.com
library(magrittr)
library(ggplot2)
theme_set(theme_minimal())

metadata <- fread("ihmp/data/ibdmdb_healthy.csv")
sra_data <- fread("ihmp/data/ibdmd_healthy_runtable.tsv")[, .(Run, `Library Name`)]
sra_data[, "External ID" := tstrsplit(`Library Name`, "_MGX")[[1]]]
metadata <- metadata[sra_data, on="External ID", nomatch=0]
metadata
man <- fread("db/data/manifest.csv")[, .(db, rank, matched_taxid, seqlength)]

foods <- fread("ihmp/data/food_abundance.csv")
foods <- man[foods, on="matched_taxid"]
foods[, "tpm" := 1.0e6 * reads / as.double(seqlength)]
foods[is.na(tpm), tpm := 0.0]

contents <- fread("ihmp/data/food_content.csv")
contents
tests <- list(
  vegetables = list(group="Vegetables (salad, tomatoes, onions, greens, carrots, peppers, green beans, etc)", ex=expression(food_group == "Vegetables")),
  fruits = list(group="Fruits (no juice) (Apples, raisins, bananas, oranges, strawberries, blueberries", ex=expression(food_group == "Fruits")),
  beans = list(group="Beans (tofu, soy, soy burgers, lentils, Mexican beans, lima beans etc)", ex=expression(food_subgroup == "Beans")),
  `white meat` = list(group="White meat (chicken, turkey, etc.)", ex=expression(food_subgroup == "Poultry")),
  shellfish = list(group="Shellfish (shrimp, lobster, scallops, etc.)", ex=expression(food_subgroup %chin% c("Mollusks", "Crustaceans"))),
  fish = list(group="Fish (fish nuggets, breaded fish, fish cakes, salmon, tuna, etc.)", ex=expression(food_subgroup == "Fishes")),
  `red meat` = list(group="Red meat (beef, hamburger, pork, lamb)", ex=expression(food_subgroup %chin% c("Swine", "Bovines", "Caprae")))
)
library(patchwork)

tables <- list()

freqs <- c(
  `No, I did not consume these products in the last 7 days` = 0,
  `Within the past 4 to 7 days` = 1/5.5,
  `Within the past 2 to 3 days` = 1/2.5,
  `Yesterday, 1 to 2 times` = 1.5,                                
  `Yesterday, 3 or more times` = 3,
  `NA` = 0
)

results <- list()
plots <- list()
merged <- metadata[foods, on=c(`Run` = "sample_id")]
for (i in 1:length(tests)) {
  name <- names(tests)[i]
  vals <- tests[[i]]
  dt <- merged[, .(abundance=sum(tpm[eval(vals$ex)]), group=.SD[[vals$group]][1], total_reads=total_reads[1]), by="External ID"]
  ns <- dt[, .N, by="group"] |> setkey(group)
  dt <- dt[group %chin% names(freqs) & ns[group, N] >= 10]
  dt[, "frequency" := freqs[group]]
  dt[, group := factor(group, levels=names(freqs)[names(freqs) %in% group])]
  dt[, "relative" := abundance]
  dt[, "item" := name]
  plots[[name]] <- ggplot(dt) + aes(x=group, y=log10(relative)) + 
    geom_jitter(width=0.3) + 
    geom_boxplot(width=0.1) + 
    stat_smooth(aes(group=1), method="lm") + 
    labs(title=name, x="consumption frequency [servings/day]", y="abundance")
  print(name)
  abmod <- lm(log10(relative + 1e-6) ~ frequency, data=dt)
  premod <- glm((abundance > 0) ~ frequency, data=dt)
  tables[[name]] <- dt
  results[[name]] <- dt[, .(n=.N, relative=mean(log10(relative*1e6)), sd=sd(log10(relative)), detected=sum(abundance > 0), item=name), by="group"]
  results[[name]][, "p_abundance" := summary(abmod)$coefficients[2, 4]]
  results[[name]][, "p_prevalence" := summary(premod)$coefficients[2, 4]]
}
[1] "vegetables"
[1] "fruits"
[1] "beans"
[1] "white meat"
[1] "shellfish"
[1] "fish"
[1] "red meat"
r <- rbindlist(results)
ta <- rbindlist(tables)
ggplot(r) +
  aes(x=group, y=relative) +
  geom_jitter(data=ta, aes(y=log10(relative*1e6)), stroke=0, size=0.5, width=0.3) +
  #geom_boxplot(outlier.color="NA") +
  geom_pointrange(aes(ymin=relative-sd, ymax=relative+sd), shape=21, fill="white") +
  geom_text(data=unique(r, by="item"), aes(x=0, y=3, label=sprintf("p[t-test]=%.2g\np[logit]=%.2g", p_abundance, p_prevalence)), size=3, vjust=1, hjust=0, nudge_y=0.5) +
  facet_wrap(~ item, ncol=2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  labs(x="", y="log-relative abundance [log₁₀RPM]") +
  guides(color=F)
Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as of ggplot2 3.3.4.

ggsave("figures/ihmp_foods.png", width=3, height=11, dpi=300)

Diet distances

Beta diversity tests

Let’s write a little function that runs the test for microbiome <-> other comparison and plots and returns the results.

Now we run it for the measures.

LS0tCnRpdGxlOiAiaUhNUCIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3J9CmxpYnJhcnkoZGF0YS50YWJsZSkKbGlicmFyeShtYWdyaXR0cikKbGlicmFyeShnZ3Bsb3QyKQp0aGVtZV9zZXQodGhlbWVfbWluaW1hbCgpKQoKbWV0YWRhdGEgPC0gZnJlYWQoImlobXAvZGF0YS9pYmRtZGJfaGVhbHRoeS5jc3YiKQpzcmFfZGF0YSA8LSBmcmVhZCgiaWhtcC9kYXRhL2liZG1kX2hlYWx0aHlfcnVudGFibGUudHN2IilbLCAuKFJ1biwgYExpYnJhcnkgTmFtZWApXQpzcmFfZGF0YVssICJFeHRlcm5hbCBJRCIgOj0gdHN0cnNwbGl0KGBMaWJyYXJ5IE5hbWVgLCAiX01HWCIpW1sxXV1dCm1ldGFkYXRhIDwtIG1ldGFkYXRhW3NyYV9kYXRhLCBvbj0iRXh0ZXJuYWwgSUQiLCBub21hdGNoPTBdCm1ldGFkYXRhCmBgYAoKYGBge3J9Cm1hbiA8LSBmcmVhZCgiZGIvZGF0YS9tYW5pZmVzdC5jc3YiKVssIC4oZGIsIHJhbmssIG1hdGNoZWRfdGF4aWQsIHNlcWxlbmd0aCldCgpmb29kcyA8LSBmcmVhZCgiaWhtcC9kYXRhL2Zvb2RfYWJ1bmRhbmNlLmNzdiIpCmZvb2RzIDwtIG1hbltmb29kcywgb249Im1hdGNoZWRfdGF4aWQiXQpmb29kc1ssICJ0cG0iIDo9IDEuMGU2ICogcmVhZHMgLyBhcy5kb3VibGUoc2VxbGVuZ3RoKV0KZm9vZHNbaXMubmEodHBtKSwgdHBtIDo9IDAuMF0KCmNvbnRlbnRzIDwtIGZyZWFkKCJpaG1wL2RhdGEvZm9vZF9jb250ZW50LmNzdiIpCmNvbnRlbnRzCmBgYAoKYGBge3J9CnRlc3RzIDwtIGxpc3QoCiAgdmVnZXRhYmxlcyA9IGxpc3QoZ3JvdXA9IlZlZ2V0YWJsZXMgKHNhbGFkLCB0b21hdG9lcywgb25pb25zLCBncmVlbnMsIGNhcnJvdHMsIHBlcHBlcnMsIGdyZWVuIGJlYW5zLCBldGMpIiwgZXg9ZXhwcmVzc2lvbihmb29kX2dyb3VwID09ICJWZWdldGFibGVzIikpLAogIGZydWl0cyA9IGxpc3QoZ3JvdXA9IkZydWl0cyAobm8ganVpY2UpIChBcHBsZXMsIHJhaXNpbnMsIGJhbmFuYXMsIG9yYW5nZXMsIHN0cmF3YmVycmllcywgYmx1ZWJlcnJpZXMiLCBleD1leHByZXNzaW9uKGZvb2RfZ3JvdXAgPT0gIkZydWl0cyIpKSwKICBiZWFucyA9IGxpc3QoZ3JvdXA9IkJlYW5zICh0b2Z1LCBzb3ksIHNveSBidXJnZXJzLCBsZW50aWxzLCBNZXhpY2FuIGJlYW5zLCBsaW1hIGJlYW5zIGV0YykiLCBleD1leHByZXNzaW9uKGZvb2Rfc3ViZ3JvdXAgPT0gIkJlYW5zIikpLAogIGB3aGl0ZSBtZWF0YCA9IGxpc3QoZ3JvdXA9IldoaXRlIG1lYXQgKGNoaWNrZW4sIHR1cmtleSwgZXRjLikiLCBleD1leHByZXNzaW9uKGZvb2Rfc3ViZ3JvdXAgPT0gIlBvdWx0cnkiKSksCiAgc2hlbGxmaXNoID0gbGlzdChncm91cD0iU2hlbGxmaXNoIChzaHJpbXAsIGxvYnN0ZXIsIHNjYWxsb3BzLCBldGMuKSIsIGV4PWV4cHJlc3Npb24oZm9vZF9zdWJncm91cCAlY2hpbiUgYygiTW9sbHVza3MiLCAiQ3J1c3RhY2VhbnMiKSkpLAogIGZpc2ggPSBsaXN0KGdyb3VwPSJGaXNoIChmaXNoIG51Z2dldHMsIGJyZWFkZWQgZmlzaCwgZmlzaCBjYWtlcywgc2FsbW9uLCB0dW5hLCBldGMuKSIsIGV4PWV4cHJlc3Npb24oZm9vZF9zdWJncm91cCA9PSAiRmlzaGVzIikpLAogIGByZWQgbWVhdGAgPSBsaXN0KGdyb3VwPSJSZWQgbWVhdCAoYmVlZiwgaGFtYnVyZ2VyLCBwb3JrLCBsYW1iKSIsIGV4PWV4cHJlc3Npb24oZm9vZF9zdWJncm91cCAlY2hpbiUgYygiU3dpbmUiLCAiQm92aW5lcyIsICJDYXByYWUiKSkpCikKYGBgCgoKCmBgYHtyLCBmaWcud2lkdGg9NCwgZmlnLmhlaWdodD0zfQpsaWJyYXJ5KHBhdGNod29yaykKCnRhYmxlcyA8LSBsaXN0KCkKCmZyZXFzIDwtIGMoCiAgYE5vLCBJIGRpZCBub3QgY29uc3VtZSB0aGVzZSBwcm9kdWN0cyBpbiB0aGUgbGFzdCA3IGRheXNgID0gMCwKICBgV2l0aGluIHRoZSBwYXN0IDQgdG8gNyBkYXlzYCA9IDEvNS41LAogIGBXaXRoaW4gdGhlIHBhc3QgMiB0byAzIGRheXNgID0gMS8yLjUsCiAgYFllc3RlcmRheSwgMSB0byAyIHRpbWVzYCA9IDEuNSwgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIAogIGBZZXN0ZXJkYXksIDMgb3IgbW9yZSB0aW1lc2AgPSAzLAogIGBOQWAgPSAwCikKCnJlc3VsdHMgPC0gbGlzdCgpCnBsb3RzIDwtIGxpc3QoKQptZXJnZWQgPC0gbWV0YWRhdGFbZm9vZHMsIG9uPWMoYFJ1bmAgPSAic2FtcGxlX2lkIildCmZvciAoaSBpbiAxOmxlbmd0aCh0ZXN0cykpIHsKICBuYW1lIDwtIG5hbWVzKHRlc3RzKVtpXQogIHZhbHMgPC0gdGVzdHNbW2ldXQogIGR0IDwtIG1lcmdlZFssIC4oYWJ1bmRhbmNlPXN1bSh0cG1bZXZhbCh2YWxzJGV4KV0pLCBncm91cD0uU0RbW3ZhbHMkZ3JvdXBdXVsxXSwgdG90YWxfcmVhZHM9dG90YWxfcmVhZHNbMV0pLCBieT0iRXh0ZXJuYWwgSUQiXQogIG5zIDwtIGR0WywgLk4sIGJ5PSJncm91cCJdIHw+IHNldGtleShncm91cCkKICBkdCA8LSBkdFtncm91cCAlY2hpbiUgbmFtZXMoZnJlcXMpICYgbnNbZ3JvdXAsIE5dID49IDEwXQogIGR0WywgImZyZXF1ZW5jeSIgOj0gZnJlcXNbZ3JvdXBdXQogIGR0WywgZ3JvdXAgOj0gZmFjdG9yKGdyb3VwLCBsZXZlbHM9bmFtZXMoZnJlcXMpW25hbWVzKGZyZXFzKSAlaW4lIGdyb3VwXSldCiAgZHRbLCAicmVsYXRpdmUiIDo9IChhYnVuZGFuY2UgKyAxKS8odG90YWxfcmVhZHMgKyAxKV0KICBkdFssICJpdGVtIiA6PSBuYW1lXQogIHBsb3RzW1tuYW1lXV0gPC0gZ2dwbG90KGR0KSArIGFlcyh4PWdyb3VwLCB5PWxvZzEwKHJlbGF0aXZlKSkgKyAKICAgIGdlb21faml0dGVyKHdpZHRoPTAuMykgKyAKICAgIGdlb21fYm94cGxvdCh3aWR0aD0wLjEpICsgCiAgICBzdGF0X3Ntb290aChhZXMoZ3JvdXA9MSksIG1ldGhvZD0ibG0iKSArIAogICAgbGFicyh0aXRsZT1uYW1lLCB4PSJjb25zdW1wdGlvbiBmcmVxdWVuY3kgW3NlcnZpbmdzL2RheV0iLCB5PSJhYnVuZGFuY2UiKQogIHByaW50KG5hbWUpCiAgYWJtb2QgPC0gbG0obG9nMTAocmVsYXRpdmUpIH4gZnJlcXVlbmN5LCBkYXRhPWR0KQogIHByZW1vZCA8LSBnbG0oKGFidW5kYW5jZSA+IDApIH4gZnJlcXVlbmN5LCBkYXRhPWR0KQogIHRhYmxlc1tbbmFtZV1dIDwtIGR0CiAgcmVzdWx0c1tbbmFtZV1dIDwtIGR0WywgLihuPS5OLCByZWxhdGl2ZT1tZWFuKGxvZzEwKHJlbGF0aXZlKjFlNikpLCBzZD1zZChsb2cxMChyZWxhdGl2ZSkpLCBkZXRlY3RlZD1zdW0oYWJ1bmRhbmNlID4gMCksIGl0ZW09bmFtZSksIGJ5PSJncm91cCJdCiAgcmVzdWx0c1tbbmFtZV1dWywgInBfYWJ1bmRhbmNlIiA6PSBzdW1tYXJ5KGFibW9kKSRjb2VmZmljaWVudHNbMiwgNF1dCiAgcmVzdWx0c1tbbmFtZV1dWywgInBfcHJldmFsZW5jZSIgOj0gc3VtbWFyeShwcmVtb2QpJGNvZWZmaWNpZW50c1syLCA0XV0KfQpgYGAKCmBgYHtyLCBmaWcud2lkdGg9MywgZmlnLmhlaWdodD0xMX0KciA8LSByYmluZGxpc3QocmVzdWx0cykKdGEgPC0gcmJpbmRsaXN0KHRhYmxlcykKZ2dwbG90KHIpICsKICBhZXMoeD1ncm91cCwgeT1yZWxhdGl2ZSkgKwogIGdlb21faml0dGVyKGRhdGE9dGEsIGFlcyh5PWxvZzEwKHJlbGF0aXZlKjFlNikpLCBzdHJva2U9MCwgc2l6ZT0wLjUsIHdpZHRoPTAuMykgKwogICNnZW9tX2JveHBsb3Qob3V0bGllci5jb2xvcj0iTkEiKSArCiAgZ2VvbV9wb2ludHJhbmdlKGFlcyh5bWluPXJlbGF0aXZlLXNkLCB5bWF4PXJlbGF0aXZlK3NkKSwgc2hhcGU9MjEsIGZpbGw9IndoaXRlIikgKwogIGdlb21fdGV4dChkYXRhPXVuaXF1ZShyLCBieT0iaXRlbSIpLCBhZXMoeD0wLCB5PTMsIGxhYmVsPXNwcmludGYoInBbdC10ZXN0XT0lLjJnXG5wW2xvZ2l0XT0lLjJnIiwgcF9hYnVuZGFuY2UsIHBfcHJldmFsZW5jZSkpLCBzaXplPTMsIHZqdXN0PTEsIGhqdXN0PTAsIG51ZGdlX3k9MC41KSArCiAgZmFjZXRfd3JhcCh+IGl0ZW0sIG5jb2w9MikgKwogIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlID0gOTAsIHZqdXN0ID0gMC41LCBoanVzdD0xKSkgKwogIGxhYnMoeD0iIiwgeT0ibG9nLXJlbGF0aXZlIGFidW5kYW5jZSBbbG9n4oKB4oKAUlBNXSIpICsKICBndWlkZXMoY29sb3I9RikKZ2dzYXZlKCJmaWd1cmVzL2lobXBfZm9vZHMucG5nIiwgd2lkdGg9MywgaGVpZ2h0PTExLCBkcGk9MzAwKQpgYGAKCgpgYGB7ciwgZmlnLndpZHRoPTUsIGZpZy5oZWlnaHQ9M30KcmFyZSA8LSBmb29kc1ssIC4obj1zdW0ocmVhZHMgPiAwKSwgcmVhZHM9c3VtKHJlYWRzKSwgdG90YWxfcmVhZHM9bWF4KHRvdGFsX3JlYWRzKSksIGJ5PSJzYW1wbGVfaWQiXQpnZ3Bsb3QocmFyZVtuPjBdLCBhZXMoeD10b3RhbF9yZWFkcywgeT1yZWFkcywgY29sb3I9dG90YWxfcmVhZHMpKSArIGdlb21fcG9pbnQoKSArIHN0YXRfc21vb3RoKG1ldGhvZD0ibG0iLCBjb2xvcj0idG9tYXRvIikgKwogIGxhYnMoeD0iI3JlYWRzIG1hcHBlZCB0byBmb29kIiwgeT0iI29ic2VydmVkIGZvb2QgaXRlbXMiLCBjb2xvcj0ibGlicmFyeSBzaXplIikgKyBzY2FsZV94X2xvZzEwKCkgKyBzY2FsZV95X2xvZzEwKCkKZ2dzYXZlKCJmaWd1cmVzL2Zvb2RfcmVhZHMucGRmIiwgd2lkdGg9NSwgaGVpZ2h0PTMpCmBgYCAKCiMjIERpZXQgZGlzdGFuY2VzCgpgYGB7cn0Kc3BlY2llcyA8LSBmcmVhZCgiaWhtcC9kYXRhL1NfY291bnRzLmNzdiIpW2dyZXBsKCJCYWN0ZXJpYXxBcmNoYWVhIiwgbGluZWFnZSldCm1pY3JvYmlvbWVfbWF0IDwtIGRjYXN0KHNwZWNpZXMsIHNhbXBsZV9pZCB+IGxpbmVhZ2UsIHZhbHVlLnZhciA9ICJuZXdfZXN0X3JlYWRzIiwgZmlsbD0wLCBmdW4uYWdncmVnYXRlID0gc3VtKQpzaWRzIDwtIG1pY3JvYmlvbWVfbWF0WywgdHN0cnNwbGl0KHNhbXBsZV9pZCwgIlNfIilbWzJdXV0KbWljcm9iaW9tZV9tYXQgPC0gYXMubWF0cml4KG1pY3JvYmlvbWVfbWF0WywgInNhbXBsZV9pZCIgOj0gTlVMTF0pCnJvd25hbWVzKG1pY3JvYmlvbWVfbWF0KSA8LSBzaWRzCmdvb2QgPC0gY29sTWVhbnMobWljcm9iaW9tZV9tYXQgPiAwKSA+IDAuNQptaWNyb2Jpb21lX21hdCA8LSBtaWNyb2Jpb21lX21hdFssIGdvb2RdCiMgbWljcm9fcmFyZSA8LSBwaHlsb3NlcShvdHVfdGFibGUobWljcm9iaW9tZV9tYXQsIHRheGFfYXJlX3Jvd3MgPSBGKSkgfD4gcmFyZWZ5X2V2ZW5fZGVwdGgoMTAwMDAwKSB8PiBvdHVfdGFibGUoKQptaWNyb2Jpb21lX3JlbGF0aXZlIDwtIG1pY3JvYmlvbWVfbWF0IC8gcm93U3VtcyhtaWNyb2Jpb21lX21hdCkKYGBgCgpgYGB7cn0KZGlldCA8LSBtZXRhZGF0YVssIG5hbWVzKG1ldGFkYXRhKVs3Mjo5Ml1bLTExXSwgd2l0aD1GXQpkaWV0IDwtIGFwcGx5KGRpZXQsIDIsIGZ1bmN0aW9uKHgpIGlmZWxzZSh4ID09ICIiLCAwLCBmcmVxc1t4XSkpCnJvd25hbWVzKGRpZXQpIDwtIG1ldGFkYXRhW1siUnVuIl1dCmRpZXQgPC0gZGlldFtyb3dTdW1zKGRpZXQpID4gMCwgXQpkaWV0IDwtIGRpZXRbLCBjb2xNZWFucyhkaWV0ID4wKSA+IDBdCmRpbShkaWV0KQpgYGAKCmBgYHtyfQpmb29kYWIgPC0gZm9vZHMgPC0gZnJlYWQoImlobXAvZGF0YS9mb29kX2FidW5kYW5jZS5jc3YiKVtyZWFkcyA+IDBdCmZvb2RhYlssICJpZCIgOj0gcGFzdGUobWF0Y2hlZF90YXhpZCwgc3BlY2llcywgd2lraXBlZGlhX2lkLCBzZXA9InwiKV0KZm9vZF9tYXQgPC0gZGNhc3QoZm9vZGFiLCBzYW1wbGVfaWQgfiBpZCwgdmFsdWUudmFyID0gInJlYWRzIiwgZmlsbD0wLCBmdW4uYWdncmVnYXRlID0gc3VtKQpzaWRzIDwtIGZvb2RfbWF0Wywgc2FtcGxlX2lkXQpmb29kX21hdCA8LSBhcy5tYXRyaXgoZm9vZF9tYXRbLCBzYW1wbGVfaWQgOj0gTlVMTF0pCnJvd25hbWVzKGZvb2RfbWF0KSA8LSBzaWRzCmZvb2RfcmVsYXRpdmUgPC0gZm9vZF9yZWxhdGl2ZVtyb3dTdW1zKGZvb2RfcmVsYXRpdmUpID4gMCwgXQpmb29kX3JlbGF0aXZlIDwtIGZvb2RfbWF0IC8gcm93U3Vtcyhmb29kX21hdCkKZm9vZF9yZWxhdGl2ZSA8LSBmb29kX3JlbGF0aXZlWywgY29sTWVhbnMoZm9vZF9yZWxhdGl2ZSA+MCkgPiAwXQpmb29kX3JlbGF0aXZlWzE6MywgMTozXQpgYGAKCiMjIEJldGEgZGl2ZXJzaXR5IHRlc3RzCgpMZXQncyB3cml0ZSBhIGxpdHRsZSBmdW5jdGlvbiB0aGF0IHJ1bnMgdGhlIHRlc3QgZm9yIG1pY3JvYmlvbWUgPC0+IG90aGVyIGNvbXBhcmlzb24gYW5kCnBsb3RzIGFuZCByZXR1cm5zIHRoZSByZXN1bHRzLgoKYGBge3J9CmxpYnJhcnkodmVnYW4pCgptaWNyb19tYW50ZWwgPC0gZnVuY3Rpb24oZmlyc3QsIG90aGVyLCBkZXNjcmlwdGlvbikgewogIHNpZHMgPC0gaW50ZXJzZWN0KHJvd25hbWVzKGZpcnN0KSwgcm93bmFtZXMob3RoZXIpKQogIG1pY3JvX2Rpc3QgPC0gdmVnZGlzdChmaXJzdFtzaWRzLCBdLCAiYnJheSIpCiAgb3RoZXJfZGlzdCA8LSB2ZWdkaXN0KG90aGVyW3NpZHMsIF0sICJicmF5IikKICAKICB0ZXN0IDwtIG1hbnRlbChtaWNyb19kaXN0LCBvdGhlcl9kaXN0LCBtZXRob2Q9InBlYXJzb24iKQogIHJlc3VsdHMgPC0gZGF0YS50YWJsZShtaWNybz1hcy5udW1lcmljKG1pY3JvX2Rpc3QpLCBvdGhlcj1hcy5udW1lcmljKG90aGVyX2Rpc3QpKQogIHBsIDwtIGdncGxvdChyZXN1bHRzKSArCiAgICBhZXMoeD1vdGhlciwgeT1taWNybykgKwogICAgZ2VvbV9wb2ludChzaXplPTEsIGFscGhhPTAuMjUsIHN0cm9rZT0wKSArCiAgICBzdGF0X3Ntb290aChtZXRob2Q9ImxtIikgKwogICAgbGFicyh4PXBhc3RlKGRlc2NyaXB0aW9uLCAiW0JyYXldIiksIHk9InNwZWNpZXMgYWJ1bmRhbmNlIGRpc3RhbmNlIFtCcmF5XSIpCiAgcHJpbnQocGwpCiAgcHJpbnQodGVzdCkKfQpgYGAKCk5vdyB3ZSBydW4gaXQgZm9yIHRoZSBtZWFzdXJlcy4KCmBgYHtyfQptaWNyb19tYW50ZWwobWljcm9iaW9tZV9yZWxhdGl2ZSwgZGlldCwgIkZGUSBkaXN0YW5jZXMiKQpgYGAKCmBgYHtyfQptaWNyb19tYW50ZWwobWljcm9iaW9tZV9yZWxhdGl2ZSwgZm9vZF9yZWxhdGl2ZSwgIk1FREkgZm9vZCBkaXN0YW5jZXMiKQpgYGAKCmBgYHtyfQptaWNyb19tYW50ZWwoZm9vZF9yZWxhdGl2ZSwgZGlldCwgIkZGUSB2cy4gTUVESSIpCmBgYAo=